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Abstract 

Efficient computation of lattice defect geometries such as point defects, dislocations, disconnections, 
grain boundaries, interfaces and free surfaces requires accurate coupling of displacements near the defect to 
the long-range elastic strain. Flexible boundary condition methods embedded a defect in infinite harmonic 
bulk through the lattice Green function. We demonstrate an efficient and accurate calculation of the lattice 
Green function from the force-constant matrix for general crystals with an arbitrary basis by extending a 
method for Bravais lattices. New terms appear due to the presence of optical modes and the possible loss of 
inversion symmetry. By separately treating poles and discontinuities in reciprocal space, numerical accuracy 
is controlled at all distances. We compute the lattice Green function for a two-dimensional model with 
broken symmetry to elucidate the role of different coupling terms. The algorithm is generally applicable 
in two and three dimensions, to crystals with arbitrary number of atoms in the unit cell, symmetry, and 
interactions. 
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I. INTRODUCTION 



Atomic-scale modeling of lattice defect geometries such as point defects, dislocations, discon- 
nections, grain boundaries, interfaces and free surfaces are key for understanding many material 
properties[l]. The anisotropic elasticity solutions for the displacement fields are accurate at dis- 
tances far from the defects (21 [31; however, the solutions frequently diverge near the defect cen- 
ter, requiring an atomistic approach to determine the defect core geometry. For many defects, 
especially dislocations, the long-range strain field is incompatible with periodic boundary condi- 
tions. To efficiently model these defects, Sinclair et al. developed a flexible boundary condition 
method based on the lattice Green function (LGF) [31 • The technique was extended for crack 
propagation[[51|6l using empirical potentials and recently applied to isolated dislocations with den- 
sity functional theory (DFT)['7|| and empirical potentials [[8l|9l. Prior to this method, dislocation 
calculations in DFT were limited to dipole ffTOl [m and quadrupole lfT2l [T3l geometries due to the 
long-range strain field of an isolated dislocation. The flexible boundary condition method couples 
the simulation cell to infinite bulk by treating an intermediate region away from the defect core 
as harmonic and relaxing these forces with a LGF. Efficient numerical calculation of the LGF for 
point defects in cubic lattices is well known [[T4l - [T6l . An automated technique for efficiently calcu- 
lating the lattice Green function with arbitrary atomic interactions was only applicable to Bravais 
lattices IfTTll . making it unsuitable for many materials with more than one atom in the crystal basis 
such as HCP (e.g. Mg and Ti). 

We extend this numerical technique to general crystals with arbitrary numbers of atoms in the 
crystal basis. The extension to multiple atom unit cells requires a separate treatment for acoustic 
and optical modes in the long wavelength limit, but continues to rely on the same input informa- 
tion (force constants) and shows similar efficiency [[TSll . Section |ll] describes harmonic response 
in a multiatom basis and the general symmetries of the dynamical matrix and lattice Green func- 



tion. Section III describes the numerical technique for efficiently calculating the LGF for general 



crystals. Section IV computes the lattice Green function for a simple doubled square lattice and 
new terms due to symmetry breaking. This symmetry breaking is manifest in both reciprocal and 
real space. The result is an efficient, automated algorithm for calculating the lattice Green func- 
tion for general crystals which can efficiently relax defect geometries; in particular, relaxation of 
dislocation geometries in Mg [[T9ll , Ti [[20ll . and as part of an interfacial LGF calculation of a twin 
boundary in Ti [[2ni . 
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II. HARMONIC RESPONSE 



Rexible boundary condition techniques allow for efficient calculation of isolated lattice defects 
with a small number of atoms by using the perfect lattice Green function (LGF) to couple the 
defect to bulk. We extend the numerical method for calculating the LGF IfTTll to work with general 
crystals with multiple atoms in the unit cell. The infinite harmonic crystal is well known from 
classical and quantum theory Il22ll23ll . For a crystal with atoms in the basis, the 3N x 3A^ force- 
constant matrix D. - determines the force on basis atom i at lattice site ^ in Cartesian 
direction a from a displacement of a basis atom j at lattice site ^' in Cartesian direction 



i2 T rtotal 



D JR-R') = 



(1) 

it=0 



duiAR)dujf^(R') 

Due to independence of differentiation order and the inversion symmetry of all Bravais lattices, 
the force-constant matrix obeys .^0^ = Dj^.^{-^). Unlike a Bravais lattice, a general crystal 
does not necessarily have inversion symmetry. Therefore, Dj^.^(-^) = Dj^.^(M) is nof guaranteed. 
The force-constant matrix obeys a sum rule, 

J]D.^jfs(R) = 0, '^i,a,/3 (2) 

due the absence of forces under a uniform translation of the crystal. Under a uniform strain, there 
is no net force on the unit cell. Hence, the force-constant matrix obeys 

J]D.^j;,iR)iR + ^i-^j)fi = 0, Va (3) 

where the 4 are the positions of the atoms within the unit cell. Since there is no torque under a 
uniform rotation, the force-constant matrix also obeys 

2 + ^' - - Diajy(R)(R + ^j)/^) = 0' V/, a,/3, y. (4) 

In the harmonic limit, the LGF G.^j^{R - R') gives displacements (uia) in response to forces (fj^), 

UiaiR) = -Yj ^iaJR - R')fdR') (5) 

where a and /3 index Cartesian directions, / and j index atoms in the crystal basis, and R and R' 
are lattice vectors. The LGF is the pseudoinverse of the force-constant matrix, 

X ^ia,iy^R-R")Gi^,jp^R" -R') = 6ij6^p6{R-R'), ^/i, j,a,/3,R,R' . (6) 



The sum rule (Eq. (|2])) guarantees that D.^j^(R) has three zero modes (uniform translation), and 
therefore D.^ j^{R) is singular. 

ni. COMPUTATION OF THE LGF 

For computational efficiency and control of numerical errors, we compute the LGF by inverting 
the dynamical matrix in reciprocal space. First, we Fourier transform the force-constant matrix to 
the dynamical matrix. We then invert the dynamical matrix using a block partitioning scheme 
by separating the dynamical matrix into acoustic and optical modes in order to isolate the poles 
and discontinuities. The inverse contains a first- and second-order poles and a discontinuity in 
reciprocal space. For numerical efficiency and stability, we perform the inverse Fourier transform 
to real space analytically for the poles and discontinuities, and numerically for the semi-continuum 
correction. Finally, to get the real space LGF, we rotate back into the crystal coordinate system to 
complete the calculation. 

Computing the LGF is more tractable in reciprocal space, where the Fourier transforms of the 
LGF (similarly for the force constant matrix and dynamical matrix) are. 



for unit cell volume V. Note that we choose the Fourier phase factor to correspond to the crystal 
vector between two atoms. We use tildes and underlines to represent matrices in reciprocal and 



The sum rule (Eq. ([2])) means that for ^ = 0, Diaj^(0) has three zero eigenvalues corresponding 

-* — -> 

to the three uniform translation modes. In addition, for small k, three modes of Dfajpik) will go 

as k ; this creates a singularity in Giaj^ik) atk = 0, corresponding to both second- and first-order 

poles. Hence, we will expand the lattice Green function around k = 0, and solve for the individual 

terms in the expansion from the expansion of the dynamical matrix around ^ = 0. 

To isolate the second-order pole in ffie elastic Green function due to translational symmetry, we 

choose a basis for atomic displacements/forces in the unit cell to separate the acoustic and optical 




BZ 



(7) 



real space, respectively. In reciprocal space, Eq. (|6]) becomes 




(8) 
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modes at ^ = 0. This involves the eigenvectors of Djaj^ik = 0), 



y, £»,-.,7.(0)4 = 



(9) 



we identify the first three jx = 1 ... 3 eigenvalues M' = as the acoustic modes, and the remaining 
3A^ - 3 as optical modes with positive eigenvalues. The acoustic eigenvectors are e^^ = d^pl^jN 
for = 1 ... 3, and the full set of eigenvectors provide an orthonormal basis. In this new basis, the 
dynamical matrix is 



<T(T',fIV 



*= z 



1 + (ik ■ (R + Xj - X,)) 



ik-iR + - f,))^ 
2! 



+ 



(10) 



where cr = A for //= 1 ... 3, and cr = O otherwise, with a similar relation between cr' and v. For 
the acoustic-acoustic projection (cr = cr' = A), the zeroth order term is zero due to the sum rule 
Eq. ©, 

^ Z = ^ Z Z ^J.'^^^^ = 0- (11) 



The remaining odd order contributions in ^to the acoustic-acoustic quadrant are zero due to inver- 
sion symmetry of a Bravais lattice (R -R). Thus, the AA term expands as 



DAA,iiv{k) = ^ 



( {k-{R + Xj- xi)f (k ■ (R + Xj - xi))^ 

+ TT- -I- 



2! 



4! 



(12) 
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Since the leading order term of the acoustic subspace of the dynamical matrix is second order in 
k, the acoustic-acoustic (AA) projection of the LGF will have a second-order pole in k. Similarly, 
the sum rule also requires the zeroth order term of the acoustic-optical (AO) and optical-acoustic 
(OA) projections of the dynamical matrix to be zero. The AO and OA projections have first-order 
poles, and the optical-optical (00) projection does not have a pole. In the rotated basis, we can 
perform a block inversion of the dynamical matrix 
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' J^AA 


Dao 
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Goo , 






Dqo ^ 



-1 



(Daa 


- Z)ao£>oo^Io) 


- (Daa - DaoDoqD^q) ^ao^oo 


-^Bo^Io 


Daa - DaoD~]jD\q) 


(Doo - Dl^Dl^^D^oY' 



(13) 



where the roman indexes A and O correspond to the projection onto the acoustic and optical 
bases respectively. This ensures that the k~^ divergence as ^ ^ is contained in the acoustic- 
acoustic quadrant of the matrix. Divergences of order ik~^ can also appear at leading order in 



the acoustic-optical and optical-acoustic quadrants for crystals without inversion symmetry. For 
crystals with inversion symmetry, only the even order terms in the acoustic-acoustic quadrant and 
the odd order terms in the acoustic-optical and optical-acoustic quadrants of the dynamical matrix 



remain (Eq. ([TO])). We write our expansion of the dynamical matrix and LGF in power series, 

DO DO 

n=0 n=-2 

where cr and cr' can be basis A or O, and D^^, and G|^^, are the nth coefficients of the power series. 
We compute the power series coefficients for the dynamical matrix by calculating the n^^ term of 



the expansion in Eq. ( fTO] ), 



- (k ■ (R + Xj - xi)Y „ ^ 

^^^.vW = Z ^ nD.,,iR)e], (15) 

As stated previously, = = 0, D% = dJ?^ = 0, G'^-^^ = G^-^^ = 0, and G^;^^ = G^^f = 0. 
In terms of these power series coefficients, the divergent and discontinuous terms of the LGF 

are 

Gif = A(^)- 

= ^(2)-l^(3)^(2)-l 
Gi°l = A(2)-lA(4)A(2)-l _^(2)-l^(3)^(2)-l^(3)^(2)-l 

(16) 

/^(-i) _ ri-m _ n(o)-in(i)t a(2)-i 

^(0) _ ^(0)t _ 75(0)-175(2)t a(2)-1 _ n(0)-iK(i) 75(0)-i75(i)t a(2)-1 , T^W-lT^dW *(2)-l * (3) * (2)-l 
^OA ~ ^AO ~ ^OO ^AO GO ^OO^OO ^AO ^ -^OO ^AO 

~(1) 



^OO - i^OO - ^AO ^AA ^Aoj 



where. 



^AA ^ ^AO^OO ^AO 

a(3) _ ?i(0)-i7i(2)t _ 7^(2) 7;(0)-i7C(i)t , 7^(1) KW-iT^d) 7^(0)-i 9^(1)1 

~ ^AO^OO ^AO ^AO^OO ^AO ^AO^OO ^OO^OO ^AO 

A(4) - _n(i) /'n(0)-in(2) n(0)-i _ nP^-iRd) n(0)-in(i) nC'^-i'l n^')"^ d?^ 

~ ^AOV oo ^oo^oo ^OO ^OO^OO ^OO^OO ^^AO ^ ^ 

, n(2) f^m-if^(2n , f^w n(0)-i7^(3)t , 7^(3) f:m-\f:iin 

^AO^OO -^AO -^AG^OO ^AG ^AO^OO ^AO 

- n(2^ fim-ifiim _ Kd) n^^^-inCD 75(0)-i 75(2)1 _ 75(4) 

^AO^GO ^00^00 ^AO ^AG^GG ^GG^GG ^AG ^AA 

For crystals with inversion symmetry, the relations for some of the terms are considerably simpli- 
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fled 



^(0) _ ^(0)t _ K(0)-175(l)t *(2)-l *(3)*(2)-l 
^OA ~ ^AO ~ ^OO ^AO^^ 

A (3) _ n(0)-i n(0)-iK(i)t 

~ ^AO^OO ^OO^OO ^AO 



(18) 



A(4) 



Kd) n(0)-i7;(3)t , 7^(3) 7^(0)-i7^(i)t _ 7^(4) 

-^AO^OO ^AO ^AO^OO -^AO ^AA 



Inverse Fourier transforming the divergent terms converges slowly. To efficiently calculate the 
LGF in real space, we integrate these divergent terms analytically, and integrate the remaining 
continuous terms numerically. To facilitate analytic integration, we expand the divergent terms 
in spherical harmonics for a 3D LGF and in a Fourier series for a 2D LGF. We apply a smooth 
spherical cut-off function /cut(^/^max) 



1 







< X < a, 
a < X < \, 

1 < X 



(19) 



where ^^ax is the radius of a sphere inscribed in the Brillouin zone, to the divergent terms so that 
they and their derivatives are zero at the Brillouin zone edge. The LGF semicontinuum correction 
is defined as the term remaining after subtracting the divergent and discontinuous terms, 

G'^A(k) = Gaa(^) - (-k-^G^^l\k) - ik-'G^;l\k) + GZ(k))U(k/k^,,) 
gS(^) = G^o(^) = Gao(^) - (-ik-'G^-^\k) + GZ(k)) U(k/k^,,) (20) 
G'ooik) = Gooik) - (GZ(k))fcui(k/k^,d 

where we treat the semicontinuum terms through numerical inversion and subtraction of the di- 



vergent and discontinuous terms from Eq. ( [T6| ). For small values of k (k < fcmax/10), numerical 
truncation error in the divergent terms dominates the calculation. Instead of direct subtraction of 
the divergent and discontinuous terms at small k, we define the leading order terms of the quad- 
rants, as 



-r 


Gaa 


-ik-^G^'^^ 

IK. Lr^Q 


-ik' 


Gqa 





R{k) = 

We then calculate the semicontinuum correction for small k as 

G' 



(21) 



G'\k) = 



[\--Hr\t.)D{ll)) 



E'\k) + ik- 



AA 







Gaa 


Gao 


Gqa 






(22) 



In 3D, a spherical harmonic expansion represents the angular variation in the Green function 
power series. The spherical harmonic coefficients are 

^max ^ 

^^^(^) = £ Z ^f.,JUh (23) 

/=0 ;ji=-/ 

where the series is truncated for / > Lmax, and L^ax is chosen such that the Im components above 
Lmax are less than 10^^^ of the largest Im component below L^ax, and Yim(J<.) are the normalized 
real spherical harmonics for k = (sin(0) cos((/i), sin(0) sin(0),cos(0)). Due to inversion symmetry 
reciprocal space, for n even, only the even / terms are nonzero, and for n odd, only the odd / terms 
are nonzero. The spherical harmonic components are found numerically by integrating over a 
uniform grid in and with Gaussian quadrature in ^ IfTTll . Rotating back with the eigenvectors e^^, 
and taking the inverse Fourier transform of the divergent and discontinuous terms gives the real 
space form lfTTll 

= Z Z ^Z)^-mY''" ,g J j, ^-^y^TTi k^'^fUm^MklR + x, - f,|), 

1=0 m=~l \\J^ + Xi-Xj\) ^-^ Jo 

(24) 

where ji(x) is the Z'*^ spherical Bessel function of the first kind, and V is the volume of the unit 
cell. The radial integrals are calculated numerically using adaptive Gauss-Kronrod integration. 
Finally, the semicontinuum term is inverse Fourier transformed using a uniform k-point mesh 
to a tolerance of 10"'. The error in the numerical inverse Fourier transform scales as A^^^'' for 
dimensionality rf[[T8||. 

For a 2D LGF, such as one used to relax an infinite, straight line defect, we expand in a Fourier 
series. The plane of the 2D LGF is normal to the threading vector f of the defect; f defines the 
periodicity of the defect. By summing along f, the Fourier transform of the LGF reduces to 2D in 
the Brillouin zone normal to L The Fourier coefficients are 

G^lik) = £ G^Z'-y-' (25) 

m=0 

where k = (cos(0), sm(9)), and the Fourier series is truncated at Mmax so that components above 
Mmax are less than 10"'^ of the largest component below M^ax- Due to inversion symmetry in k 
space, the coefficients for even (odd) n are only nonzero for even (odd) m. Rotating back with the 
eigenvectors e^j^, and taking the inverse Fourier transform of the divergent and discontinuous terms 
of the 2D LGF gives the real space form lfTTl 

^Sjp^^) = E '^S?;^™^ "'^(-1)"'^ V k'^"fcAk/km.Mk\^ + Xi- fyl), (26) 

m=0 ^0 



where Jmix) is the m* Bessel function of the first kind, and A = V/\i\ is the area of the 2D unit cell. 
We then evaluate these Bessel integrals numerically using adaptive Gauss-Kronrod integration. 
Finally, the semicontinuum term is inverse Fourier transformed using a uniform k-point mesh to a 
tolerance of 10"'. 

In summary, we: (1) Determine the acoustic/optical basis by calculating the eigenvectors e^^ 
of the dynamical matrix atk = 0. (2) Fourier transform the force-constant matrix and rotate the 



dynamical matrix into the acoustic/optical basis (Eq. ([TO])). (3) Determine the power series coeffi- 
cients of the dynamical matrix ^^(k) (Eq. ([T5])) on an angular grid for analytic block inversion 



of the divergent and discontinuous terms of the LGF G''^^, ^^(k) (Eq. 16 and Eq. 17|). (4) Deter- 



mine the 2D Fourier coefficients G^^, (Eq. ([25])) or 3D spherical harmonic coefficients G^^]., 



(Eq. (|23])). (5) Calculate the semi-continuum term G^^ik) on a regular grid in reciprocal space by 
interpolating the Fourier or spherical harmonic expansion onto the grid and subtracting from the 



block inverse of the dynamical matrix (Eq. (20) and Eq. (22)). (6) Inverse Fourier transform the 



divergent and discontinuity terms into real space analytically (Eq. ( |26] ) or (|24])), numerically in- 
verse Fourier transform the semi-continuum terms. (7) Add all contributions and rotate back to the 
original atomic basis to get the LGF in real space. This automated algorithm directly calculates 
the lattice Green function from the force constant matrix with controllable numerical errors for a 
general crystal with any number of atoms in the crystal basis. The numerical error scales with the 
number of k-points A^^. as N^^'', where d is the dimension of the system lfTSll . In the next section, 
we determine leading order corrections to the lattice Green function in terms of force-constant 
symmetry breaking. 



IV. SQUARE LATTICE MODEL 
A. Single atom unit cell 

Fig. [T] is a schematic of the doubled, elastically isotropic square lattice with ffist- and second- 
nearest-neighbor radial springs of spring constants 1 and 1/2 (^ = 77 = 0) that illustrates the 
additional terms added by a multiple-atom basis. The single atom lattice with lattice vectors 
di = (a, 0) and 02 = (0, a) has elastic constants 

3 11 

Cii = C22 = -, C12 = C21 = -, Cgg = -• (27) 
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FIG. 1 . Schematic of the two atom basis square lattice nearest neighbor interaction model with the spring 
constant perturbations ^ and rj labeled. The unit cell is shaded in gray. Two different asymmetries are 
considered: 77 describes the asymmetry between the white-white and black-black atom interactions (along 
y), and ^ describes the asymmetry between the white-black unit cell spring and the white -black neighboring 
cell interactions (along x). The diagonal second neighbor "bond-bending" springs of strength 1 /2 are shown 
in red; these springs stabilize the lattice and produce isotropic elastic response when rj = ^ = 0. 

This material is elastically isotropic (ICge = Cn - Cu) with radial interactions for Cauchy sym- 
metry (C12 = Cee); or 2D Poisson ratio of 1/3, Young's modulus of 4/3, and shear modulus of 1/2. 
The Fourier transform of the LGF for the single atom case is shown in Fig.|2j The divergent and 
discontinuous terms in the LGF for the single atom square lattice are 



3k^a 



^ 72 



2 - cos 26 - sin 26 
- sin 26 2 + cos 26 

7 -4cos26'-hcos46' 

7 -h 4 cos 20 -h cos 40 



(28) 



B. Doubled unit cell 

Fig. [T] shows the doubled square lattice with lattice vectors = (2a, 0) and a2 = (0, a) and 
crystal basis: Xb = (0,0) and x„ = (a,0) (for "black" and "white" atoms). We consider the same 
interactions as in the single atom case. The Fourier transform of the two atom LGF is shown in 
Fig.|3j The acoustic-acoustic quadrant, which corresponds to the summed interactions of the black 
and white atoms, is the same as the LGF for the single atom case except for the halved Brillouin 
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FIG. 2. Contour plot of the Fourier transform of the lattice Green function for an isotropic square lattice 
= T] = 0) over the first Brillouin zone (BZ). The LGF is Hermitian, thus the bottom-left is the Hermitian 
conjugate of the upper-right. All components show a second-order pole a.tk = 0; c.f. Fig. HI 



zone. The divergent and discontinuous terms in the LGF are 



PT(-2) ^ ^ 



2 - cos 26 
- sin 26 



- sin 26 
2 + cos 26 



r(-^) _ ri-m _ a 



r(0) - J_ 

AA 72 



7 -4cos 26* -h cos 40 





7 -I- 4 cos 26* -I- cos 46* 



^(0) _ 1 



1 

3 



(29) 



with all other quadrants zero, and the same elastic constants as the single atom case. 

Fig. |4] shows polar plots of the divergent and discontinuous LGF terms. Inversion symmetry 
requires the acoustic-optical quadrant to be purely imaginary; since this is a doubled single atom 
system, the acoustic-optical quadrants are zero. The acoustic-optical quadrants correspond to the 
response of the internal degrees of freedom of the system to elastic strain. The doubled system 
behaves just as the single atom system; thus, there is no internal relaxation. The second-order 
pole and the discontinuity correction in the acoustic-acoustic quadrant are the same as the single 
atom case. This is because both cases describe the long range elastic behavior summed over all 
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FIG. 3. Contour plot of the Fourier transform of the lattice Green function for a doubled isotropic square 
lattice = ?7 = 0) over the first Brillouin zone (BZ) in the acoustic(A)/optical(0) rotated basis. The BZ 
is cut in half along the direction due to the doubling of the lattice along the x direction. Doubling the 
cell also results in the appearance of optical quadrants in the LGF. The LGF is Hermitian, thus the bottom- 
left triangle is the Hermitian conjugate of the upper-right triangle. The acoustic-acoustic quadrant, which 
corresponds to the collective motion of all the atoms, corresponds to the single atom unit cell LGF, and has 
a second-order pole at ^ = 0; c.f. Fig.|4j 

atoms. The leading order optical-optical constants correspond to the inverse of the optical phonon 
frequencies at the Y point and are not discontinuous. Since the cell is doubled along x, the doubled 
system is stiffer along the xx mode than the yy mode because the xx mode involves both bond- 
stretching and -bending springs, but the yy mode only involves bond-bending springs. 
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FIG. 4. The directional dependence of the leading order divergent and discontinuous terms of the lattice 
Green function for a doubled square lattice. The elastic Green function corresponds to the 1 jk^ term. The 
optical-optical quadrant is not discontinuous because the leading order corresponds to the inverse of the 
optical modes at ^ = 0. The magnitude along k in the polar plots is the multiplier for the LGF for direction 
k; dashed lines correspond to negative multipliers. 



C. Breaking single atom symmetry 



We introduce small perturbations which break the single-atom symmetry of the doubled lattice 
to produce changes in the LGF. We break the symmetry of the black/white atoms by changing the 
black-black spring constant to (1 +7]) and the white- white spring constant to (1 - 77) as in Fig.[T} 
This has no effect on long range elastic behavior, and so the poles of the LGF are unmodified 
by this perturbation. To leading order in 77, the interaction adds a discontinuity correction in the 
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acoustic-optical mode, 

/ 

pr(0) _ ^(0)T _ J]_ 

^OA ^AO J 2 





2 sin 26*- sin 40 -3 + 4 cos 20 + cos 40 



+ 0(77'). (30) 



The acoustic-acoustic discontinuity term is modified at second order in 77. Fig. [5] shows polar 
plots of these perturbations. The acoustic-optical quadrants of the LGF correspond to the response 
of internal modes to collective motion. However, breaking the black/white atom symmetry alone 
is not enough to generate internal relaxation in response to long range elastic waves as there is 
no adjustment to the l/k^ or i/k poles of the LGF. The leading order adjustment occurs in the 
discontinuity correction which only depends upon the direction of k. Thus, there is an angular- 
dependent change due to the different stiff"nesses of the black-black and white-white interactions. 

To break symmetry and introduce internal relaxation, we set the in-cell black-white spring 
constant to (1 -1- ^) and the black- white spring constant to neighboring cells to (1 - ^) as shown 
in Fig. [T] This spring constant change causes a change in the long range elastic behavior of the 
model. The Cn elastic constant becomes, to second order, 3/2-^^/12 due to internal relaxation of 
the atoms in the unit cell in response to strain, while C22 remains 3/2. Linear order in ^ introduces 
the acoustic-optical i/k pole 

„\ 

+ 0(e). (31) 



'(-1) _ pr(-i)t _ 2^' cos 9 



-2 COS 20 sin 20 




At second order in ^, there are changes in the acoustic-acoustic second-order pole and disconti- 
nuity correction, and the optical-optical discontinuity correction. Fig. [5] shows polar plots of the 
first-order perturbations. By breaking the in-cell symmetry, an imaginary term now appears in the 
acoustic-optical quadrants of the LGF. This i/k pole comes from internal relaxation of atoms in 
the unit cell due to a long- wavelength elastic wave. Breaking the internal symmetry causes the 
black and white atoms to respond diff"erently to a long range strain and shift from their strained 
simple square lattice sites in response. This term is important for describing the internal response 
of a multi-atom basis crystal and does not appear in the single atom Bravais lattice case. 



D. The simple square lattice in real space 



Fig. [6] shows the real space LGF for the single and doubled unit cells in crystal space by taking 
the inverse Fourier transform. Without interaction perturbations, both are identical to 10"^ for a 
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FIG. 5. Leading order divergent term corrections of the lattice Green function for a doubled square lattice 
as a function of the scaled relative spring constants rj (black-black/white-white coupling) and f (black- 
white/white-black coupling). The magnitude along ^ in the polar plots is the multiplier for the LGF for 
direction k; dashed lines correspond to negative multipUers. 

40 X 80 k-point mesh for the doubled cell and a 80 x 80 k-point mesh for the single cell. Since both 
have the same interactions, both methods yield the same harmonic response in a different basis. 
The LGF plot shows the linear relationship between force and displacement for atoms separated 
by ^+ ^ - Xj. The elastic LGF is logarithmic in \R-\-Xi- Xj\ and dominated by the inverse transform 
of the k'^ pole at long range. 

To first order, neither the rj or ^ perturbations modify the pole of the LGF. Therefore, at 
long range in real space, the LGF is still dominated by the original isotropic elastic behavior. Con- 
sidering the inverse transforms of the perturbed interactions, the rj perturbation results in opposite 
yy response for the black/black and white/white interactions as expected since it is a change in 
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FIG. 6. LGF in real space for an isotropic square lattice = 77 = 0). The LGF is only defined at crystal 
sites. Red and blue coloring correspond to negative and positive values of the LGF at those crystal sites. 
The two atom LGF is identical to the single atom LGF at the crystal sites to 10^^ for a 80 x 80 (40 x 80) k 
point mesh for the single (double) atom cell. 

the yy spring strengths. The ^ perturbation results in opposite xx response for the black/white and 
white/black interactions. The effects of ^ are seen in the ik'^ pole, and only visible at short to 
intermediate range with a decay of 1 /R, while rj first appears in fc" component of the LGF, and is 
only seen at short range since it decays as 1 /R^. Figure [v] shows the first-order real-space response 
to ^ and 77 perturbations. 



V. CONCLUSION 

The direct, automated algorithm can efficiently and accurately calculate the lattice Green func- 
tion for general crystals with more than one atom in the unit cell basis and arbitrary interactions. 
Additional terms describing the response of internal degrees of freedom of the system correspond- 
ing to optical modes appear in this formalism that did not appear in a treatment for a Bravais 
lattice. Including the additional optical terms of the lattice Green function extendes the previ- 
ous automated calculation of the LGF for long-range interactions IfTTll to general crystals. This 
technique efficiently calculates defect structures in HCP metals such as Mg [|T9l . Ti EOl . as well 
as semiconductors and intermetallics using flexible boundary condition methods[7]. In particu- 
lar, reducing the number of atoms required for accurate calculation of isolated dislocation core 
geometries provides efficient use of density functional theory. 
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FIG. 7. Change to the black atom LGF due to the ^ (black-white/white-black coupling) and rj (black- 
black/white-white coupling) perturbations. The LGF is only defined at crystal sites, with red and blue 
coloring corresponding to negative and positive values. The change in LGF for white atoms is opposite to 
the change for black atoms. 
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